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Method for Measuring Structural Thickness from 
Low-Resolution Digital Images 

CROSS REFERENCE TO RELATED APPLICATIONS 

[0001] This application claims priority to US Provisional Application No. 60/43 1,129 on 
December 5, 2002, the content of which is herein incorporated in its entirety. 

GOVERNMENT INTEREST 

[0002] This invention was supported in part by Grant Nos. ROl 41443 and R21 471 12 from 
the National Institutes of Health. Accordingly, the Government may have certain rights in this 
invention. 

FIELD OF THE INVENTION 

[0003] The invention deals with the extraction of object features from digital images 
acquired at low resolution, specifically, the measurement of structural thickness distribution 
along an object. Targeted appUcations comprise, but are not limited to, the measurement of 
trabecular bone thickness in magnetic resonance or computed tomography images. 

BACKGROUND OF THE INVENTION 

[0004] Effective tools for object shape analysis are important and useftil in many imaging 
applications including medical ones. One such popular and widely used tool is the distance 
transform (DT) (Rosenfeld et al. Pattern Recog. 1:33-61 (1968); Danielsson, Comput, Graphics 
Image Process. 14:227-248 (1980); Borgefors, Comput. Vision Graphics Image Process. 
27:321-345 (1984); Borgefors, Comput. Vision Image Understanding 64:368-376 (1996)) of an 
object. For a hard (binary) object, DT is a process that assigns a value at each location within 
the object that is simply the shortest distance between that location and the complement of the 
object. 

[0005] However, until very recently this notion of hard DT have not been applicable on 
fiizzy objects in a meaningfiil way (Kaufinann, Introduction to the Theory of Fuzzy Subsets , 
Vol. 1, Academic Press, New York, 1975; Rosenfeld, Inform. Control 40:76-87 (1979); Bezdek 
and S. K. Pal, Fuzzy Models for Pattern Recognition . IEEE Press, New York, 1992). The 
notion of DT for fuzzy objects, called fuzzy distance transform (FDT), becomes more important 
in many imaging applications because it is often necessary to deal with situations with data 
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inaccuracies, graded object compositions, or limited image resolution (on the order of an 
object's structure size). 

[0006] The notion of fuzzy distance is formulated by first defining the length of a path on a 
fiizzy subset and then finding the infimum of the lengths of all paths between two points. For 
the fiizzy distance between two points x and j/, the space n is defined as the "infimum" 
(Weisstein, CRC Concise Encyclopedia of Mathematics , Chapman & Hall/CRC, Boca Raton, 
FL, 1999) of the lengths of all paths between them. The length of a path tt in a fiizzy subset of 
the w-dimensional continuous space 9?" is defined as the integral of fiizzy membership values 
along 71. Generally, there are an infinite number of paths between any two points in a fuzzy 
subset and it is often not possible to find the shortest path, if it exists. Thus, the fiizzy distance 
between two points is defined as the infimum of the lengths of all paths between them. It is 
demonstrated that, unlike in hard convex sets, the shortest path (when it exists) between two 
points in a fiizzy convex subset is not necessarily a straight-line segment. For any positive 
number 9 < 1, the 0 -support of a fiizzy subset is the set of all points in 91" with membership 
values greater than or equal to 9 . It is shown that, for any fiizzy subset, for any nonzero 6 < 1, 
fiizzy distance is a metric for the interior of its 9-support. 

[0007] The FDT is thus defined as a process on a fiizzy subset that assigns to a point its 
fiizzy distance fi*om the complement of the support. The theoretical firamework of FDT in 
continuous space is extended to digital cubic spaces and it is shown that for any fiizzy digital 
object, fiizzy distance is a metric for the support of the object. In general, FDT is useful, for 
example, in feature extraction (Fu et aL, IEEE Trans, Comput. 25:1336-1346 (1976)), local 
thickness or scale computation (Pizer et al, Comput. Vision Image Understanding 69:55-71 
(1998); Saha et al, Comput Vision Image Understanding 77:145-174 (2000)), skeletonization 
(Srihari et al,, in Proceedings of International Conference on Cybernetics and Society^ Denver, 
Colorado, pp. 44-49 (1979); Tsao etal, Comput. Graphics Image Process, 17:315-331 (1981); 
Saha et al.. Pattern Recog. 30:1939-1955 (1997)), and morphological (Serra, Image Analysis 
and Mathematical Morphology , Academic Press, San Diego, 1982) and shape-based object 
analyses (Borgefors, "Applications of distance transformations," in Aspects of Visual Form 
Processing (Arcelli et al, eds.), pp. 83-108, World Scientific, Singapore, 1994). In particular, 
FDT may be usefiil in fault detection in integrated circuit chips or in computer motherboard 
circuits, analysis of the dynamics of a hurricane, etc. FDT will be usefiil in many medical 
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imaging applications, such as computation of local thickness of trabecular bone or vessels, or 
morphology-based separation of anatomic structures having similar intensities, e.g., artery-vein 
separation. 

[0008] Trabecular or cancellous bone - the type of bone that dominates in the vertebrae and 
at locations near the joints of long bones (metaphysis and epiphysis) - consists of a network of 
plates and struts. Bone atrophy as it occurs in osteoporosis, leads to either homogeneous or 
heterogeneous thinning of the trabecular elements. Besides changes in network connectivity 
(and thus of the topology of the network) the thickness of the trabeculae most critically 
determines the mechanical competence and thus resistance to fracture of trabecular bone. 
Accurate measurement of trabecular thickness is, therefore, of significant interest, for example, 
to assess the effectiveness of anabolic (bone forming) agents of patients with osteoporosis. 
[0009] The classical approach toward measuring trabecular thickness has been based on 
histomorphometry of transiliac bone biopsies (Chavassieux et aL, in Osteoporosis, 2 (Marcus et 
a/., eds.) New York: Academic Press, pp. 501-509 (2001)). Typically, the perimeter of the 
trabeculae is measured in stained sections, and thickness is approximated as the bone area 
divided by one half of the perimeter (Parfitt et al, J. Clin, Invest 72:1396-409 (1983)). The 
emergence of imaging technologies, such as micro computed tomography (|i-CT) (Ruegsegger 
et aL, Calcified Tissue International 58:24-29 (1996)) enables reconstruction of three- 
dimensional images calling for more elaborate techniques for measuring structural thickness. 
[0010] One model-independent approach involves inscribing spheres into the structure 
(Hildebrand et aL, J. Microscopy 185:67-75 (1997)) in such a manner that trabecular thickness 
at any location is computed as the diameter of the largest inscribed sphere containing that 
location. The implementation issues are solved using distance transform and the distance ridge, 
which provides the set of the center points of largest inscribed spheres. This approach is well 
suited for high-resolution images that can easily be segmented, but it is bound to fail when 
significant partial volume blurring is present. The latter is the case in images acquired in the 
limited spatial resolution regime of in vivo |x-MRI and ji-CT that are beginning to supplant bone 
biopsy-based methods for structural analysis of trabecular bone (see, e.g., WehrH et aL, Topics 
in Magnetic Resonance Imaging 13:335-356 (2002). 

[001 1 ] However, the fuzzy nature of these images, caused by partial volume blurring, 
virtually precludes binarization. Accordingly, until the present invention a long felt need has 

-3- P-2944 

520631 1 



Attorney Docket 22253-73413 



remained in this art for a method that obviates segmentation and that can effectively deal with 
images acquired at a voxel size comparable to the typical trabecular bone thickness. It is a goal 
therefore, to better understand the fuzzy distance transform (in both continuous and digital 
spaces), to study its properties, to present a dynamic programming-based algorithm to compute 
FDT for fiizzy digital objects, and to demonstrate practical applications of the FDT methods. 

SUMMARY OF THE INVENTION 

[0012] The present invention first provides the theory underlying the algorithm in the digital 
cubic space and the algorithm itself, and then further provides utilitarian demonstrations of the 
method for extracting structural thickness from fiizzy distance transform images. The 
robustness underlying the present FDT-based methods is shown in terms of CT images of 
human trabecular bone in low-resolution regimes, re-sampled at progressively coarser 
resolution, and after application of rotation and addition of noise (to simulate the in vivo 
situation). Reproducibility of embodiments of present methods is demonstrated on the basis of 
images from |i-CT volume data by comparing thickness histograms from individual sections 
from the same and from different subjects, as well as with MRI volume data sets from human 
volunteers imaged consecutively in different sessions. Using in vivo ji-MR images from prior 
animal studies, it is shown that short-term exposure to a bone atrophy-inducing agent results in 
trabecular thinning that can be quantified with a method of the present invention. 
[0013] The theoretical framework of FDT is developed in the n-dimensional continuous 
space n. Consider a path n over a fuzzy subset of n. Unlike the case of hard sets, the 
membership values of the points through which n passes need to be considered to determine its 
length. Thus, in an embodiment it is illustrated that, unlike the case of hard convex sets, the 
shortest path (when it exists) between two points in a fuzzy convex subset is not necessarily a 
straight-line. For any positive number 0 < 1, the 0-support of a fiizzy subset is the set of all 
points in 9i" with membership values greater than or equal to 9. 

[0014] In a further embodiment, it is shown that for any fuzzy subset, for any nonzero 6 < 1 , 
fiizzy distance is a metric for the interior of its 9-support. 

[0015] It is also shown that, for any smooth fiizzy subset, fiizzy distance is a metric for the 
interior of its 9-support (referred to as "support"), FDT is defined as a process on a fiizzy subset 
that assigns to a point its fiizzy distance from the complement of the support. The theoretical 
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framework of FDT in the continuous space is extended to cubic digital spaces, and in a further 
embodiment it is shown that for any fuzzy digital object, fuzzy distance is a metric for the 
support of the object. 

[0016] In another embodiment of the invention, a dynamic programming-based algorithm is 
presented for computing FDT of a fuzzy digital object. It is shown that the algorithm terminates 
in a finite number of steps and when it does so, it correctly produces the FDT image. 
[0017] Further, exemplary applications for fuzzy distance transforms in medical imaging are 
presented, including the quantification of the thickness of vascular structures and of trabecular 
bone over a wide range of resolution regimes from |j,-CT to in vivo |x-MRI data, performance of 
the method under object rotation, sensitivity to noise, reproducibiUty of thickness values, and 
efficacy. A resolution-dependent correction of average thickness has been shown to be effective 
in largely eliminating the systematic error resulting from failure to sample FDT values at the 
true skeleton locations. Further, FDT-derived thickness was found to be remarkably robust to 
noise, and the values were stable for SNR > 5. 

[0018] In addition, in an embodiment of the method has been shown to be reproducible in 
repeat MRI scans conducted in vivo demonstrating its applicability for longitudinal studies for 
the assessment of drug efficacy. Presented are several potential, exemplary applications of 
fuzzy distance transform in medical imaging, including quantification of blood vessels and 
trabecular bone thickness in the regime of limited special resolution where these objects become 
fuzzy. 

[0019] Additional objects, advantages and novel features of the invention will be set forth in 
part in the description, examples and figures which follow, all of which are intended to be for 
illustrative purposes only, and not intended in any way to limit the invention, and in part will 
become apparent to those skilled in the art on examination of the following, or may be learned 
by practice of the invention. 

BRIEF DESCRIPTION OF THE DRAWINGS 

[0020] The foregoing summary, as well as the following detailed description of the 
invention, will be better imderstood when read in conjunction with the appended drawings. It 
should be imderstood, however, that the invention is not limited to the precise arrangements and 
instrumentalities shown. 
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[0021] FIG. 1 depicts an example demonstrating that the shortest paths (when they exist) 
between two points in a convex fiizzy subset are not necessarily a straight-Hne segment. The 
support (S) of the fuzzy disk is the union of different shaded regions. The points within the 
dark gray region have a sufficiently high membership value in S while the points within the 
light gray region have a sufficiently low membership value. The shortest path between the two 
, points X and is contained within the light gray region and is not a straight-line segment. 
[0022] FIG. 2 depicts an example demonstrating why a raster scan based approach fails to 
compute FDT for fuzzy digital objects in a fixed number of scans. The support 0(0) of the 
fuzzy object O is the union of different shaded regions. The dark (light) gray region resembles 
the set of points with sufficiently high (respectively, low) membership values in O. The path n\ 
(712) is the shortest path in O from p\ (respectively, /?2) to the complement of 0(0). 
Computation of the length of n\ (712), /.e., the FDT value at p\ (respectively, 772) in O, needs three 
(respectively, five) raster scans. 

[0023] FIGs. 3(a)-3(c) depict the application of FDT on a fuzzy object. FIG. 3(a) shows a 
2D slice from a 3D CTA image of a patient's head. FIG. 3(b) shows a fuzzy object representing 
the bone structures in FIG. 3(a). The fuzzy object was constructed from the image of FIG. 3(a) 
by first thresholding the bone regions, and then by blurring and subsequently adding noise to it. 
FIG. 3(c) shows a FDT image of the fiizzy object in FIG. 3(b). As shown in this figure, the 
ridges of FDT values follow the medial axis of the fuzzy object. 

[0024] FIGs. 4(a)-4(c) depict application of FDT-based thickness computation to an arterial 
tree. FIG. 4(a) shows a MEP rendering of a 3D subvolume taken from a 3D CTA image of a 
patient's head (after removing bones), showing a portion of the carotid arterial tree. FIG. 4(b) 
shows a MIP projection of the fuzzily segmented arterial tree. FIG. 4(c) shows a MIP 
projection of the FDT image of the 3D image shown in FIG. 4(b). Mean and standard deviation 
of thickness values computed along the curve skeleton of the arterial tree mask are 2.74 mm and 
1.8 mm, respectively. 

[0025] FIGs. 5a and 5b depict high-resolution in vivo 3D MR images by the FDT-based 
thickness computation method of the human distal radius showing trabecular bone structure. 
FIG. 5a shows a 2D shoe taken from the raw 3D MR image. The central highlighted disk 
represents the cross section with the cylindrical ROI used for computing thickness of the bone 
trabeculae. FIG. 5(b) shows a 3D projection of the BVF image computed over the ROL Mean 
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and standard deviation of the thickness values of trabecular bone over the selected ROI are 1 02 
//m and 42 /^m, respectively. 

[0026] FIGs. 6(a)-6(f) demonstrate preliminary validation of FDT-based thickness method 
using ^'CT data of human distal radius trabecular bone by studying thickness distributions over 
four regions, two each from the same trabecular bone sample. FIGs. 6(a), 6(b) show a fuzzy 
shell rendering of the 3D BVF image of trabecular bone networks for the two different 
specimens. FIGs. 6(c), 6(d) show two representative slices of the raw ju-CT image, one from 
each of the two samples shown in FIGs. 6(a) and 6(b). FIG. 6(e) shows two FDT-based 
thickness distributions, one from each of the two different slices (330 jum apart) chosen from the 
sample shown in FIG. 6(a). FIG. 6(f) is the same as FIG. 6(e), but for the sample in FIG. 6(b). 
[0027] FIG. 7 is a graphic illustration to demonstrate how digitization introduces a 
resolution-dependent systematic negative error in FDT-based average thickness computation. 
The shaded object region represents a cross section of trabecular bone over a small region. The 
ideal skeleton of these trabeculae is indicated by the dotted lines, and the true thickness would 
be obtained if FDT values could be sampled along the actual skeleton. However, since a digital 
skeleton is a set of grid points close to the ideal skeleton, thickness distribution is computed by 
sampling FDT values at these grid points only, which introduces a random negative error (small 
white bars). 

[0028] FIGs. 8(a)-8(d) depict image sets. FIGs. 8(a) and 8(b) depict representative 
transverse ju-CT images of human distal radius trabecular bone specimens; while FIGs. 8(c) and 
8(d) depict 3D surface-rendered BVF images. In the remaining Figures and discussion, the 
bone sample displayed in FIG. 8(c) is referred to as "sample A," and the cross-section in FIG. 
8(a) is "slice A." Those displayed in FIGs. 8(b) and 8(d) are referred to as "sample B" and 
"slice B," respectively. 

[0029] FIGs. 9(a)-9(f) depict images at comparative resolution values. FIGs. 9(a) and 9(b) 
show slices A and B, as above, from the BVF images of samples A and B at 22 jum isotropic 
voxel size, respectively. FIGs. 9(c) and 9(d) show sUces from the BVF images at 88 jum 
isotropic voxel size from samples A and B matching to slices A and B, respectively. FIGs. 9(e) 
and 9(f) show the same as FIGs. 9(c) and 9(d), respectively, but at 176 //m isotropic voxel size. 
[0030] FIG. 10 graphically depicts plots of average thickness values at different resolutions 
with and without corrections for samples A and B. Note that without correction, the computed 
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thickness values underestimate 'true' thickness with the magnitude of the error increasing 
approximately linearly with voxel resolution. 

[003 1 ] FIGs. 1 1 (a)- 1 1 (c) depict images of a slice from sample B for three different rotation 
angles around its z-axis. FIG. 1 1(a) = 35°; FIG. 1 1(b) = 125° and FIG. 1 1(c) = 245°. 
[0032] FIG. 12 graphically depicts plots of average thickness values at different in plane 
rotations around the z-axis (the axis along the slice direction). Rotation-thickness plots show 
periodic peaks at integral multiples of 45°. 

[0033] FIGs. 13(a)- 13(c) depict slice A after adding zero mean, multiplicative random noise 
at different SNR. FIG. 13(a) (SNR=1). FIG. 13(b) (SNR-5(b)). FIG 13(c) (SNR= 10). 
[0034] FIG. 14 graphically depicts the performance of the FDT method under varying 
signal-to-noise ratios. 

[0035] FIG. 15 graphically depicts two pairs of FDT-based trabecular bone thickness 
distributions (without resolution dependent corrections). Each is for a pair of spatially close 
slice locations from samples A and B. 

[0036] FIGs. 1 6(a)- 1 6(d) demonstrates the reproducibility of the FDT-derived trabecular 
bone thickness determinations for in vivo situations. FIG. 16(a) shows a slice from an in vivo \X' 
MR image of a human distal radius acquired at 137 x 137 x 350 lam^ voxel size. The actual 
region for which thickness was computed was manually outUned; the boundary of this region is 
shown in white. FIG. 16(b) shows the region of interest zoomed from the scan in FIG. 16(a). 
FIG. 16(c) shows a BVF image derived from image within the VOI shown in FIG. 9(a). FIG. 
16(d) shows a FDT image computed from the BVF-image in FIG. 16(c). 
[0037] FIGs. 17(a)-17(f) depict visually matched cross-sectional images (distal femoral 
epiphysis) of a dexamethasone treated rabbit studied at three time points: FIG. 17(a) = baseline, 
FIG. 17(b) = 4-weeks, and FIG. 17(c) = 8-weeks after onset of treatment. FIGs. 17(d) and 17(e) 
show computed BVF images corresponding to those in FIGs. 17(a)- 17(c). It is noted that no 
corticosteroid-induced trabecular bone thinning is visually apparent. 

[0038] FIG. 18 graphically depicts the effect of corticosteroid exposure on rabbit trabecular 
bone. "Sham" represents the sham-operated (untreated) group, while "treated" represents the 
corticosteroid-treated group. Plots represent mean average thickness values (after 
resolution-dependent correction) for three time points described in FIG. 17, i.e,, baseline, 4- 
weeks, and 8-weeks. Vertical bars represent standard deviations. 
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DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS OF THE INVENTION 

[0039] The present invention provides theory and algorithms for fuzzy distance transforms 
(FDT), representing an extension of the concept of distance transforms for hard objects to the 
more common fuzzy objects present in imaging. The section entitled Theory describes the 
theory and properties of FDT in both continuous and digital cubic spaces, and further describes 
the metric property of fuzzy distance. The section entitled Algorithm describes a dynamic 
programming-based algorithm for computing FDT of digital objects. The termination of the 
algorithm in a finite number of steps, as well as its correctness are thus established. The results 
of applications of the FDT method are presented in the Examples. 

[0040] It has been shown that fuzzy distance is a metric for the interior of the support of a 
fuzzy object. A saUent feature of fuzzy objects is the property that the shortest path between 
two points is no longer a straight-line. It has also been shown that the raster scan algorithm 
commonly used for distance transforms of hard objects fails in the case of fuzzy objects. A 
dynamic programming-based algorithm is presented in the invention for computing FDT of 
fuzzy digital objects. It has been shown that the algorithm terminates in a finite number of steps 
and when it does so, it correctly computes FDT. Finally, several potential applications for fuzzy 
distance transforms in medical imaging are presented, including the quantification of the 
thickness of vascular structures and of trabecular bone to voxel size, object rotations, and signal- 
to-noise ratios (SNR). In addition, reproducibility of the methods was assessed in specimen and 
in vivo studies in humans, and changes in trabecular thickness were evaluated in response to 
specific drug treatment. 

THEORY 

[0041] In this section, the theory is presented and the properties of fuzzy distance transform 
of fuzzy subsets are defined either in the continuous space or in a digital space. First described 
is the case for the continuous space, which will then guide the formulation in a digital space. 
7. FDT in Continuous Space 

[0042] Let denote the ^-dimensional continuous space. A "fuzzy subset" (see, e.g,, 
Kaufinann, 1975) S of is defined as a set of pairs {{x, |is(jc)) | x e 9?"} where |is : 9?" [0, 
1] is the "membership function" of S in 91". For any value 0 e [0, 1], "0 -support" of S, 
denoted by 0e(S), is the hard subset {x\xe y{"}3nd [is (x) > 0 }of 9?". In other words, the 0 - 
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support of S is the set of all points in 91" with membership values greater than or equal to 0. 0- 
support will often be referred to as "support" and be denoted by 0(5). A fiizzy subset with a 
bounded support is called "bounded." 

[0043] The following notions on fuzzy subsets are used herein, (see, Rosenfeld, Fuzzy Sets 
Systems 13:241-246 (1984); Bogomobiy, Fuzzy Sets Systems 23:257-269 (1987) for details). A 
fuzzy subset S is a "ring" if [i^ = |I(r) , where r = || .x - xo || for some xq e 9?" and p[ : 91 -» [0, 
1] is a membership function. S is said to be "convex" if, for every three coUinear points y, 
and z in 91^ such that;; lies between jc and z, |as (y) ^ min[|is (jc), \is (z)]. A convex ring is called 
a "fuzzy disk." A fuzzy subset is called "smooth" if it is differentiable at every locations e 91". 
[0044] Let S be a (hard) subset of 91". S is used to denote its complement, while "Interior 
(S)" denotes its interior, which is the largest open set contained in S. The distance transform 
(DT) of S may be represented as an image {(x, Ds (x)) \ x e 9?"} on 9i" where Ds is the DT 
value at x that is defined as follows. 



where, inf gives the infimum of a set of positive numbers and || • || is the Euclidean norm. 
[0045] In the above equation, S should be a proper subset of 9?". In digital images, bounded 
objects are always used so that S is always nonempty. In the subsequent discussions, the fuzzy 
distance transform is considered of only bounded subsets. 

[0046] Here, FDT of a fiizzy subset S is defined in 91". Similar to ordinary DT, FDT is an 
image on 9?". The FDT image is denoted by a set of pairs {{x, 0,six)) \ x e 9?'^}, where 0.s(x) is 
the FDT value at x which is defined in the following way. 

[0047] A "path" tc in 91" firom a point x e 91" to another point (not necessarily distinct) y e 
91" is a continuous function n : [0, 1] -> 91", such that 71(0) = x and 7r(l) =y. The "length" of a 
path n in S, denoted by Il^n), is the value of the following integration 



[0048] Following the above equation, Tlsin) is the integral of membership values (in S) 
along 71. Note that, if the inverse path is defined as n(t) = 7i:(l-0 of tt, it can be shown that 




(Equation 1) 




(Equation 2) 
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ns(7i) = ns(K'). A question arises at this point on weighting the path length using membership 
values. Although other membership-based weights, e.g., absolute derivatives, may be useful in 
some applications, membership values themselves are used herein as weights as it is intended 
that the distance between two points is defined as the minimimi material that has to be traversed 
to proceed from one point to the other. Thus, when a path passes through a low density (low 
membership) region, its length increases slowly and the portion of the path in the complement 
of the support of S contributes no length. This approach is useful to measure regional object 
depth, object thickness distribution, etc. 

[0049] In some applications, it may be useful to consider lengths of connected paths (paths 
entirely contained in 0(S)) only. This can be conceived by incorporating a little change in the 
definition of path length as follows: 



[0050] It can be shown that both lis and TVs lead to identical FDT images and Theorems 1 
and 2 (see below), which originally proved fl^ are also true for 11' s lis is used herein to define 
fuzzy distance. 

[005 1] The set of all paths firom a point jc e 91" to another point >^ e 91'' is denoted by P(jc, 
y). It may be noted that P(x, y) contains infinitely many paths. The shortest path firom x e^i^ 
toy e 91'' in S is a path 71;^^ ^ P(x, y) such that 11^.(71;^^) < Ils{n), Vtc e P(x, y). It is worth 
mentioning that there may or may not be only one shortest path between two points in a fuzzy 
subset, and when it does exist, it may not be unique. 

[0052] The following is a 2D example where no shortest path exists between two points. 
Consider a fuzzy subset of 9?^ with its support as a disk and let us pick any two points jc and y 
near its center. Within the support of the fuzzy subset, there are only two membership values - 
a high membership value for the points on the straight-line segment joining x and y, and a low 




(Equation 3) 



where 




(Equation 4) 
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membership value elsewhere. It is not difficult to see that the straight-line segment is not 
the shortest path between jc and y and for any other path n between the two points it is possible 
to find a shorter path by further straightening the path. Although the existence of one shortest 
path between two points in a fuzzy subset is not guaranteed, the infimum of path lengths always 
exists and is unique. 

[0053] Let (x, y) denote a subset of positive real numbers defined as 4s(^, y) = {05(71) | n 
e P(x, y)} ; Le, , ^s(x, y) is the set of all possible path lengths in S between x and y. The "fuzzy 
distance" from jc e 9?" to >^ € 91" in S, denoted as (o^x.y)^ is the infimum oft,(x,y); i.e.. 



[0054] Fuzzy distance satisfies metric properties as stated in the following theorems (see 
proofs, Saha et al. Computer Vision and Image Understanding 86:171-190 (2002)). 
[0055] THEOREM 1 . For any fuzzy subset S of for any nonzero positive number 9 < 1 
fuzzy distance cos is a metric for Interior(&Q(S)). 

[0056] THEOREM 2. For any smooth fuzzy subset S ofSR^, fuzzy distance cos is a metric 
forInterioT{®Q(S)). 

[0057] Interestingly, the shortest paths (when they exist) in a fiizzy subset S between two 
points X, >^ e 9? are not necessarily a straight-line segment, even when S is convex. For 
example, consider the fuzzy disk S illustrated in FIG. 1 . In the figure, the support 0S of the 
fuzzy disk is the shaded region. In this example, within the support 0S, there are two possible 
membership values - points with high membership value are shown as dark gray and those with 
low membership value (within the outer annular region) are shown as Hght gray. Consider the 
two points X and y as illustrated in the figure. Assuming the high membership value is 
sufficiently close to one and the low one close to zero, the shortest path between x, y should be 
contained within the light gray region, and therefore is not a straight-line segment. 
[0058] Maintaining the same spirit as in the distance transform of hard sets, the FDT value 
Qs{x) at a point jc e 91'* is equal to the fuzzy distance between x and the closest (with respect to 
C05) point in 0(S) . In other words, the value of value Qs (x) is defined as follows 



ai^x,y)^inf^^x,y) 



(Equation 5) 



Q5(x) = inf {o)5(^,;;)|:^e 0(S)}. 



(Equation 6) 
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2. FDT in Digital Cubic Spaces 

[0059] Fuzzy distance transformation are next described in digital cubic spaces guided by 
the formulation of the same in the continuous space as discussed in the previous section. To 
make the description more precise, the concepts are redefined in digital spaces corresponding to 
the previous definitions in the continuous space. In order to avoid xmnecessary heaping of 
notations, the same notations are used in this section as were used in Section 1, although their 
exact meaning in a digital space may be different from that of the continuous space. 
[0060] A "digital space" D is an ordered pair (G, a), where G is the underlying digital grid 
and a is a binary relation on G that indicates the adjacency relationship between every two 
points in G. In general, a "digital grid" is a set of points in 9?", such that the interpoint distances 
are bounded below, and within a bounded distance of any point in 91" there is at least one point 
in the grid. However, most imaging systems acquire images in cubic grids, and these grids are 
simple to describe as well as to understand. Therefore, the following description is made in 
terms of cubic grids. An ^-dimensional cubic grid may be constructed by dividing 9?'^ into 
hypercuboids with n orthogonal families, each of equally spaced parallel hyperplanes. The set 
of the centers of these hypercuboids generates a cubic grid, and it is not difficult to see that, 
under a proper coordinate system, these points represent the points in where Z is the set of all 
integers. Z" is used to represent an n-dimensional cubic grid. The notion of adjacency is useful 
to define a path in a digital space and the boundary of a digital object. 
[0061] Although researchers have used fuzzy adjacencies (Rosenfeld, Pattern Recog. Lett, 
2:311-317 (1991); Udupae^ a/., Graphical Models Image Process. 58:246-261 (1996)), the 
interpretation of fuzziness of adjacencies in the context of a path is not clear. Therefore, the 
description is confined to hard adjacencies only. In other words, a:Z"xZ"->{01}. Two 
points /?, ^ G Z" are called "adjacent" if, and only if a (p, q^) = 1. In the rest of this paper, a can 
be standard 4- or 8-adjacency in 2D, 6-, 18-, or 26-adjacency in 3D, and their higher 
dimensional analogs. 

[0062] Two adjacent points are often referred to as "neighboring" points to each other. 

[0063] A "digital object" "O" is a fiizzy subset {(p, //o(p)) |p g 2"}, where //q: Z"-> [0, 1] 

specifies the membership value at each point in the object. It should be noted that, in general, 

an imaging system (Cho et al,. Foundations of Medical Imaging , Wiley, New York, 1993) 

acquires images containing information of a target object, often along with other co-objects. 
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However, following the spirit of the problem addressed here, digital objects are defined instead 
of digital images (Rosenfeld et aL, in Digital Picture Processing , 2nd ed., Vol. 2, Academic 
Press, New York, 1982). In other words, the proposed framework assumes that the target object 
is already (fuzzily) segmented from an acquired digital image using an appropriate segmentation 
method (Rosenfeld et a/., 1982, Pal et al. Pattern Recog. 26:1277-1294 (1993); Udupa et al, 
(eds.), 3D Imaging in Medicine , CRC Press, Boca Raton, FL, 1991; Sonka et al. Image 
Processing, Analysis, and Machine Vision , 2nd ed., PWS Publishing, Brooks/Cole, Pacific 
Grove, CA, 1999; Kong et al. Topology Appl 46:219-262 (1992)). This process of extracting 
an object from an image, commonly referred to as image segmentation, has been studied 
extensively for decades. Since the subject itself is an open area of research, rather than 
becoming involved with that issue, the present invention is described assuming the above 
definition of a digital object. 

[0064] The support 0(O) of a digital object O is the set of all points in 2" each having a 
nonzero object membership value, /.e., 0(0) - {p\p s and /io(p) 5^ 0}. A "path" tt in a set S 
of points from a point p e Sto another (not necessarily distinct) point q e S is a sequence of 
points p=puP2y '>*,Pm^g, such that pi e S for all 1 < / < w and pj is adjacent to pj+i for all 1 < 
j<m. The length of the path is m. 

[0065] Although only hard adjacency relations are considered to define a path, the 
Euclidean distance between two adjacent points should be used in a meaningful way in defining 
its length. A set of points Swill be called "path connected" if, and only if, for every two points, 
p.qeS, there is a path in S fi*om p to q, P(p, q) will denote the set of all paths in from a 
point /? € 2" to another point q e For the purpose of defining the length of a path, the notion 
of a link and its length is used. A "link" is a path consisting of two points. The length of a link 
{p, q) in O may be defined in different ways, e.g., (1) max{//o(p) fioiq)} ^ \\p-q% (2) V2{jio{p) + 
Mo{q)) X \\p-^\l etc. Notably, in both examples, the length of a link has two components - one 
coming from the membership values at p and at q, and the other firom the distance between the 
two points. In the disclosed embodiment of the invention, the second choice is followed for the 
length of a link. Theoretical requirements for a valid length function for links are (1) the length 
of a link (p,p) is "0," (2) the length of the links <p, q) and {q,p) |p, ^ e are equal, and (3) the 
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length of a link <p, where p € 0(0) and p^q. is greater than zero. The "length" Tloiu) of a 
path 7t = <p , . . , j^m = ^) is the sum of the lengths of all links on the path, ue.. 



Ttp^ e Pip^q) is one "shortest path" from /?€Z" to qreZ" in O if rioCTTp,^) ^ OoCtc), Vti € 
P(p,q). Unlike the case of the continuous case, one shortest path always exists between two 
points in a bounded digital fiizzy object. 

[0066] PROPOSITION 1 . For any digital space D = (Z", a) for any fuzzy object O on D 
with a bounded support, for any two points p, q g Z", there exists one shortest path from p to q. 
[0067] Proof Consider a path 7c = <pi,p2, pi, pii, pi2, pii, pi, pi+i, pm> with a 
repetition at pi, forming a loop (a path with the common starting and finishing point). 
Obviously, the length of the path <pi, p2, . . ., Pi, Pi+i, . . Pm) obtained by removing the repetition, 
is less than or equal to that of n. Therefore, to prove this proposition it is shown that among the 
paths without any repetition, there is one with the smallest length. Since the support of O is 
bounded, there are only finitely many links with nonzero values, and therefore, there are finitely 
many length values for the paths without any repetition. Hence, there is the minimum length for 
the paths between p and q and a path with the minimum length is one shortest path. 
[0068] Although, the existence of one shortest path is guaranteed in fiizzy digital objects, as 
in the continuous space, there may be multiple shortest paths between two points in a fixzzy 
digital object between two points. In the rest of this disclosure, fiizzy digital objects with 
bounded supports are considered. The fiizzy distance, from p e Z" to q e Z" in O, denoted as 
®o(P> q)^ is the length of any shortest path in O from p to q. Therefore, 



[0069] Moreover, it is shown that, similar to the case of fiizzy subsets in the continuous 

space, fiizzy distance in a digital object defined over a digital cubic space is a metric. 

[0070] THEOREM 3. For any digital cubic space D = (Z", a) for any digital object O, fuzzy 

distance coo is a metric for 0(0). 

[0071] For proof of Theorem 3, see Saha et aL, 2002. 



(Equation 7) 



(Equation 8) 
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[0072] Following the framework of FDT for fuzzy subsets in the continuous space, the FDT 
value at a point g 2" in a fuzzy object O over a digital space, denoted by Qo(p), is equal to the 
fuzzy distance between p and the nearest point in 0(0) . In other words, the value of Q^(p), is 
defined as follows: 

iP) = ^®o iP^ <l) (Equation 9) 

ALGORITHMS 

[0073] In this section, an algorithm embodiment is presented for computing the FDT of 
digital objects. Also, a proof will be provided that the algorithm terminates in a finite number of 
steps, and when it does so it outputs the FDT of the digital object. Assuming a uniform 
neighborhood, we can think that, with respect to a point /? € Z", all its adjacent neighbors are 
ranked. A vector 5, called a "resolution vector,'' is used whose fth element gives the continuous 
distance between a point and its fth adjacent neighbor. This uses the (3'^-l)-adjacency relation 
{i.e., 8-adjacency in 2D and 26-adjacency in 3D). Therefore, 6 is a (3''-l)-dimensional vector. 
The information about the resolution vector may be directly obtained from the application 
imaging system. Let iV*(p) denote the set of adjacent neighbors of a point /? e Z" excluding p 
itself 

[0074] As demonstrated by Borgefors in 1996, a raster scan approach effectively computes 
regular DT for binary images using only two scans. A basic reason behind this is the fact that in 
a binary image the shortest path from a point to the background (the complement of the binary 
object) is always a straight-line segment (in a digital sense). However, this is not true for fuzzy 
digital objects (see FIG. 8). This makes a raster scan based approach inappropriate in 
computing FDT, as illustrated by example in FIG. 8. Similar to the example shown in FIG. 7, in 
FIG. 8, dark gray indicates high membership in the object O and light gray indicates low 
membership. Consider two points p\ and p2 as shown in the figure. Let ti\ be a shortest path 
from p\ to 0(0 ) , and let 712 be the same, but from pi. It is not difficult to see that three raster 
scans are needed to compute the length of 711 in O (/.e., the FDT value at p\), while it takes five 
scans to do so for m (/.e., the FDT value at pi)^ Therefore, the number of raster scans needed to 
compute the FDT of a fuzzy object is dependent on the shape of the object as demonstrated in 
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FIG, 2. The iterative algorithm (Rosenfeld et al,y 1982) for computing the distance transform of 

a hard digital object can be used to compute FDT of fiizzy objects. However, such approaches 

will be inefficient in this application as the shortest path fi-om a point to the background may be 

quite complicated, and for each iteration the method needs to scan the entire image. 

[0075] In the following section, a dynamic progranmiing-based algorithm is presented to 

compute FDT of a fuzzy digital object. 

[0076] ALGORITHM. Compute_FDT. 

Input: O = (Z^, fi ^, a, and 5, as defined above. 

Auxiliary Data Structures: modified fiizzy object O = (Z", //q), an image (Z", Q) 
representing FDT of O, and a queue Q of points. 
Output: an image (Z", Q) representing FDT of O. 
L input 0 = (Z", fio)\ 

2. for all p € 0(0) , set Q(p) =0; 

3. for all p e 0(0), set Q(p) 

/* MAX is a large value */ 

4. for all p G 0(0), such that N^(p) n 0(0) is non-empty 

5. push p into Q; 

6. while Q is not empty do 

7. remove a point p from Q; 

8. find distmm = min^6Ar*(p) [Q(^) + drank(p^) X y2(Mo(p) + juo(g))]'y 

/* rankipy q) gives the rank of ^ in the neighborhood of /?*/ 

9. find c/wtmin < ^{p) then 

10. set Q(p) = rfr5tniin; 

1 1 . push all points q € N*(p) n 0(0) into Q\ 

12. output the FDT image O 

[0077] In the following two propositions, it is proven that the embodied algorithm 
compw^e_FZ)r terminates in a finite number of steps, and when it does so, it produces the FDT 
image. 

[0078] PROPOSITION 3.1. For any fuzzy object O = (Z", //o) over any digital cubic space^ 
(Z", a), the algorithm compute JFDT terminates in a finite number of steps, 
[0079] Proof The algorithm is iterative, and at each iteration in the while-do loop, it 
removes exactly one point fi"om queue Q. Also, only the points of 0(0) visit Q, Since the 
number of elements in 0(0) is always finite, the algorithm compute_FDT fails to terminate in a 
finite number of steps only when some point p € 0(O) is modified infinitely many times. 
Following Steps 9 and 10, the value at p strictly decreases after every modification. Moreover, 
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following Steps 8, 9, and 10, during every modification,/? is assigned the value of the length of 
some path in O from /? to a point in0(O) . First, it is shown that p is always set to OoCti) for 
some 71 without any repetition or loop. If this is not true, without loss of generality, let p be the 
first point to break the rule. Therefore, the path ti is of the form qu qi^ But, this is 

possible only when p has previously been assigned the length of the path v!-q\, qi, p. 
Following Equation 7, it is obvious that no(7r') ^ no(7i). But this contradicts Step 9. Now, 
since 0(0) is finite, there are only finitely many length values for repetition-free paths from 
0(0) to p. Hence, any point p could be modified only finitely many times. Hence 
compw/e_FZ)7' terminates in a finite number of steps. 

[0080] To prove the correctness of the algorithm the following notation is needed. Let 
r 0(0) (0 denote the set of all paths of less than or equal to / points from p to some q e 0(0 ) . 

[0081] PROPOSITION 2. For any fuzzy object O = (Z" , ja^ ) over any digital cubic space 
(Z", a), when the algorithm compute JFDT terminates, its output equals the FDT image 

(z^Qo). 

[0082] Proof. Initially, the algorithm starts with the points of 0(0) (see Steps l-A), Thus, 
Steps 8-10 guarantee that whenever a point p is set to c/wtmin (in Step 10), there is always a path 
from some q e 0(0 )Xop with length disirtm^ Therefore, for any /? e Z", Cl{p) > Qq (/?) . 
[0083] It may be observed from the algorithm that it never increases the value of any /? e Z". 
It is now proven that, for any /? € Z", Qq {p) > Q(/?) by using the method of induction. 

Obviously, for anyp € Z", T^-^^i)) is nonempty implies that p e 0(0) , and following 

Equation 7, the length of any element of T^^^{\) is always "0." Thus, in Step 2, it is 

guaranteed that the algorithm sets every point /; € Z" to the length of the shortest path of 

r ^ (1) in O . Now assume that the algorithm sets every point p € 0(0 ) to the length of the 

shortest path of T ^(i -1) , for some / > 1, in O . It is shown that the algorithm sets every 

point /? € Z" to the length of the shortest path of T^-^^ii) in O in Steps 7-11. Let 

= {q = Pi, P2, Pi = p) be the shortest in O among all paths of r^(z) from a 
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point q e 0(0) to p. Following the hypothesis of the induction method, at a certain iteration 
the algorithm sets the point pi.\ to Q( ) = min Flo (tt) < Uq «/?, , /?2 > • • • > Pi-i )) Step 

10, and Step 1 1 ensures that, in the same iteration,/? is pushed into Q; and upon the removal of p 
from Q, it is assigned (see Steps 7 to 10) a value less than or equal to Q(p) = Q(p,-i) + ^rankipi^pi ,) 

X Vz (MoiPi) + (MoipiA)) < no (7Cp) (following Equation 7); note that S.^^^^^^.^. =|| ||). 

Therefore, for any p e O.o(p) > Qoip). Hence, by the results in the previous paragraph, for 
any/> G Z", Q.o(p) = Qo{p). 

FDT'BASED THICKNESS COMPUTATION 

[0084] Thickness is a useful parameter in analyzing object shape and morphology, and as 
shown above, FDT provides the depth mapping at each point within the support of a fuzzy 
object O. The principle of FDT-based thickness computation is to sample these depth values 
along the medial axis (or medial surface) of target O; thus, providing the regional thickness 
distribution over the object. "Skeletonization" is a widely used technique to generate the medial 
siu-face representation of a digital object (Srihari et al, 1979; Tsao et al, 1981; Saha et al, 
1997). This is accomplished by reducing a 3D object into a union of surfaces and curves (some 
investigators have further reduced it to a union of curves only (Pothuaud et aL, J. Microsc. 
199:149-161 (2000))), thereby preserving the topology and shape of the original object. 
[0085] Although, topology preservation of digital objects has a generally accepted 
definition, the same is not true for shape, and different researchers often use different 
definitions. The inventors' previously described method (Saha et al.. Pattern Recognition 
30: 1939-1955 (1997)) has been used in the present work for skeletonization. The method 
iteratively erodes an object along its border while preserving its topology and retaining the 
surfaces and curves necessary to describe the shape of the object. Each iteration consists of 
three complete scans of the entire image space, and in each scan there are different types of 
border points (referred to as "s-open," "e-open," and "v-open" points). The simple point 
characterization, (see also, Saha et al.y IEEE Transactions on Pattern Analysis and Machine 
Intelligence 16:1028-1032 (1994)) was used for topology preservation, i.e„ only simple points 
were deleted during erosion. Thus, efficient computational approaches were implemented. 
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[0086] To preserve object shape, edge points of surfaces and end points of curves were 
saved during erosion. Thus, both types of points are referred to as "shape points." The 
skeletonization process takes 45-60 seconds for a 3D trabecular bone image set of 256 x 256 x 
60 image size on a 450 MHz Pentium III PC running under Linux OS. The robustness of the 
skeletonization method under rotation and noise was previously demonstrated (Saha et aL, 
1997). After calculating the skeleton of the support of trabecular bone, the distribution of 
regional thickness over O was computed in accordance with the present invention by sampling 
the FDT values along the skeleton. Following the same theory, thickness is computed along the 
skeleton of an object. 

[0087] While it may be desirable to generate the skeleton of an object directly from its fuzzy 
representation by taking into account the bone volume fraction value at each location, the 
significant research and developmental demands needed are beyond the scope of the present 
invention. Therefore, instead, the standard skeleton on the support 0(0) of object O is used to 
obtain its medial axis representation and compute its thickness. It should be noted that no 
exphcit threshold is required to define the support 0(0), and therefore, the resulting skeleton is 
not threshold sensitive. 

[0088] Let denote the skeleton Sfe(0(O) of 0(0). The mean thickness of an object denoted 
by t(0), is computed as the average of twice (2x) the FDT values along the skeleton of iS'A(0(O), 
/.e., 

E 2Q^(/7) 

^ ^ ^ p^si^ (Equation 10) 

\Sk{e(p))\ 

where | S^(0(O)|, yields the number of points in the skeleton. 

[0089] As demonstrated above, the FDT approach deals effectively with partial volume 
effects. However, digitization introduces a resolution-dependent, systematic error in the FDT- 
based thickness computation, as illustrated in FIG. 1. Consequently, the present invention 
further provides an approximate solution that corrects for these digitization errors. As described 
in the previous paragraph, the true thickness can be obtained by sampling the FDT values along 
the ideal medial axis of a structure. Obviously, any deviation from this ideal medial axis can 
only reduce the estimated value of thickness. See, FIG. 1, wherein the ideal skeleton of the 
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shaded object region is shown by the dotted Unes, and the true thickness would be obtained if 
FDT values could be sampled along these dotted lines. However, in most digital image 
processing approaches, object representations, as well as related computations, are confined to 
grid points. Following this basic limitation, a digital skeleton is a set of grid points close to the 
ideal skeleton. Thus, FDT values are sampled only at grid points close to the ideal skeleton, 
thereby introducing a random negative error in computed thickness as illustrated in the figure. 
[0090] For a given voxel size, these errors have an upper bound, which is the maximum of 
the shortest distance between the ideal skeleton and a grid point. Obviously, this upper bound 
equals half of the voxel's diagonal. For an isotropic structure, assuming that the tme skeleton 
can pass through any point within the voxel with a uniform probability, the expected average 
error should be the value of the following integral: 

errors— \ \ \Jx^+y^-^z^dxdydz (Equation 11) 

where "p" is the voxel resolution. In the case of highly anisotropic bone, such as in the case of 
the distal radius, where the trabeculae are predominantly oriented along the z direction, the third 
dimension can be ignored. Using numerical approaches (see, e.g.. Press et aL, Numerical 
Recipies: The Art of Scientific Computing. Cambridge, London: Cambridge University Press, 
1986), the value of the above integration in the xy plane approximates to 

error = 03Sp 

Finally, an in-plane resolution-dependent correction was added to the average FDT-derived 
thickness values to obtain final thickness, as follows: 

Z 2ao(p) 

^j^^^ peskom ^^.o.38p (Equation 12) 

\Sk{®(0))\ ^ y J 

After execution of all preprocessing steps (BVF computation), the complete FDT based 
thickness computation method takes on the order of 1-1.5 minutes for typical 3D image data 
sets. 

[0091] To validate the results of the present FDT-based method, its performance at different 
resolutions was compared with a previously reported 3D thickness computation method (the 
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algorithm of Hildebrand et al, J. Microscopy 185:67-75 (1997)). The method estimates 
regional thickness of trabecular bone using binary distance transform and local maxima ("center 
points of all nonredundant spheres") in thresholded images. For all resolutions, a threshold of 
50% BVF was used. In addition, at the highest resolution (22 \xm) thresholds were varied 
between 25 and 75 % BVF. 

[0092] To evaluate methods of the invention and to demonstrate the potential in longitudinal 
in vivo studies, experiments were conducted to evaluate, e,g.^ the sensitivity of the thickness 
computation method to voxel size, object rotations, and signal-to-noise ratios (SNR), 
respectively, as well as reproducibility, as assessed in specimen and in vivo studies in humans, 
including specific studies to confirm the value of the method for assessing the changes in 
trabecular thickness in response to corticosteroid drug treatment. The following illustrate a few 
examples of applications of FDT methods of the present invention, and describe a method of 
computing local thickness using FDT. However, these examples are not meant to limit the 
scope of the invention, and altematives may be utilized. 

EXAMPLES 

[0093] The following experiments were all conducted in 3D, except the reproducibility 
study based on |i-CT data. For all experiments, the following steps were applied: (1) generation 
of a BVF map from the raw images, (2) FDT computation, (3) skeletonization, and (4) thickness 
computation. These steps for each experiment are described in the order in which they were 
executed. 

Example 1 - Apphcation of FDT on a fuzzv image . 

[0094] The fuzzy object was constructed from a 2D slice of a 3D computed tomographic 
angiography (CTA) image of a patient's head. The original slice is shown in FIG. 3(a). The in- 
plane resolution of the image is 0.25 x 0.25 mm^. A fuzzy object was generated as follows. The 
bone regions were segmented from the slice image through interactive thresholding. The fuzzy 
object representing the bone region was constructed by (1) blurring (to simulate partial volume 
effects) the thresholded image over a Gaussian kemel of radius 5 pixels, and (2) subsequently 
adding a correlated, zero-mean Gaussian noise with standard deviation equal to 10% of the pixel 
bone fraction value. The final fuzzy object is shown in FIG. 3(b). 
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[0095] The motivation behind the construction of the fuzzy object in this way is simply to 
demonstrate the resuhs on a relatively realistic image. FIG. 3c shows the FDT image as derived 
from the fiizzy object of 3(b). In FIG. 3(c), the intensity values are shown to be proportional to 
its FDT value. As visually apparent in this figure, the ridges of FDT values follow the medial 
axis of the fuzzy object suggesting that there is potential for applications of FDT in computing 
skeletons of fuzzy objects. 

Example 2 - Computation of thickness . 
[0096] To compute thickness in this experiment, the skeletonization method described in 
Saha et al, 1997, was used. Let Sk{@{0)) denote the skeleton of 0(0). At any point p € 
Sk{@{0)\ the thickness value was computed as twice the largest FDT value in the neighborhood 
of p. 

[0097] A first example of FDT-based thickness computation is illustrated in FIG. 4. FIG. 
4(a) shows a maximum intensity projection (MIP) rendition of a portion of an arterial tree in a 
3D CTA image of a subject's brain vasculature. The size of the image domain is 74 x 217 x 40 
mm^ and the voxel size is 0.32 x 0.32 x 1.25 mm^. The image has been rendered after removing 
bone using a recently developed method (Saha et al, in Proc. SPIE: Med, Imag., San Diego, 
CA, Vol. 4322, pp. 1264-1272 (2001)) from the inventors' laboratory. The mask for the arterial 
tree was segmented from the rest of the tissue using scale-based fuzzy connectedness (Saha et 
aLy 2000). A membership value //^ (p) at each location within the segmented artery mask was 
computed as: A. 



where /is the image intensity function, and are mean and standard deviations of intensity 



Ca as its mean and standard deviation parameters. A zero membership value is assigned at each 
location outside the artery mask. A MIP rendition of the membership image of the arterial tree 
is presented in FIG. 4(b). The FDT image was computed for the 3D membership image of the 
arterial tree and the MIP rendition of the FDT image is illustrated in FIG. 4(c). In order to 
compute thickness, a curve skeleton of the arterial tree was computed using the method in Saha 
et aLy 2000). The mean and the standard deviation of thickness (vessel diameter) values along 




(Equation 13) 



values within the arterial mask and G, 



is an unnormalized Gaussian function with and 



-23- 



P-2944 



520631 1 



Attorney Docket 22253-73413 



the skeleton were 2.74 mm and 1 .8 mm, respectively. Arterial plaques can be detected by 
identifying sudden reduction of vessel diameter while tracing along the curve skeleton of an 
arterial tree. 

[0098] A second example of thickness computation, illustrated in FIG. 5, is a high- 
resolution in vivo 3D magnetic resonance (MR) image of the human distal radius showing the 
trabecular bone network. The in vivo images were acquired on a 1.5T GE clinical scanner. The 
image size was 512 x 256 x 32 //m^ and the voxel size was 137 x 137 x 350 ^m^. A slice from 
the raw image is shown in FIG. 5(a). A cylindrical region of interest (ROI) was chosen for 
analysis. The central highlighted disk in FIG. 5(a) is the cross section of the ROI with that slice. 
[0099] The image within the ROI was preprocessed by deshading and noise reduction using 
a histogram deconvolution method (Hwang et aL, Internat J. Imaging Systems Tech. 10:186- 
198 (1999)) to produce a bone volume fraction (BVF) map. The spatial resolution was 
enhanced to 68 x68 x88 //m^ using subvoxel classification (Wu et al. Magnetic Res. Med. 
31:302-308 (1994)). A 3D projection of the final BVF image is presented in FIG. 5(b). The 
FDT image was computed from the resolution-enhanced BVF image. Thickness values were 
computed using a surface skeleton (Saha et aL, 2000) of the bone mask. The mean and the 
standard deviation of thickness values along the skeleton were 102 and 42 pim, respectively, 
in good agreement with the known thickness of human bone trabecular, which is on the order 
100-150 //m (Aaron et al. Clinical Orthopaedics Related Res. 215:260-271 (1987)). 
[0100] A further demonstration of the effectiveness of the proposed FDT-based thickness 
com-putation method involving high-resolution micro-computed tomography (//-CT) images is 
presented in FIG. 6. The basic idea was to evaluate the FDT-based thickness distributions over 
four regions, two each from the same trabecular bone sample. Two fx-CY images of different 
metaphyseal samples of human distal radius were acquired on a SANCO Medical //-CT 20 
scanner at 22 //m isotropic resolution (although other devices and at an appropriate level of 
resolution could be used). 

[0101] Each of the two //-CT images was processed as follows. The raw gray-scale image 
was binarized to yield bone masks. A BVF map at each location within the bone mask was 
computed using an equation similar to Equation 13. The BVF images of the two samples are 
illustrated in FIGs. 6(a) and 6(b) using fuzzy shell rendering (Udupa et al, IEEE Comput, 
Graphics AppL 13(6):58-67 (1993)), supported by the 3DVIEWNIX system (Udupa et al, Proc. 
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SPIE 2164:58-73 (1994)). Two rq)resentative slices of the raw //-CT images, one from each 
sample, are shown in FIGs. 6(c) and 6(d). A pair of slices, separated by 330 jum (/.e., 15 slices 
apart), was selected from each of the two samples and 2D FDT-based thickness was computed 
separately for each BVF slice image. FIG. 6(e) shows two distributions of the thickness values, 
one for each of the two different slices chosen from the sample shown in FIG. 6(a). The mean 
and standard deviations of thickness values in one sUce were 119 jum and 59 //m, respectively, 
and the corresponding values for the other slice in the same sample were 116 //m and 62 //m, 
respectively. FIG. 6(f) presents thickness distributions for the two slices chosen from the other 
sample shown in FIG. 6(b). The mean and standard deviations of thickness values in one slice 
are 103 //m and 57 //m, respectively, versus values for the other slice in the same sample of 101 
jum and 54 jum, respectively. As demonstrated in FIGs. 6(e) and 6(f), the thickness distributions 
for the slices from the same sample are quite similar while they are significantly different for the 
pair of slices from the other sample. 
Example 3 - BVF Imaging. 

[0102] For this and the experiments that follow, BVF object generation utilized CT images 
that had been acquired previously on a Scanco Medical |i-CT 20 scanner at 22 )im isotropic 
resolution from samples cored from cadaveric bone of the distal radial metaphysis of two 
donors. As above, shell rendering was by 3DVIEWNIX (Udupa et aU 1994). The samples 
represent cylindrical cores of 9 mm height and 9 mm diameter with the cylinder axis parallel to 
the direction of the radius. Because of their high spatial resolution exceeding trabecular 
thickness by a factor of 5, the data represent the "gold standard." 

[0103] Two cross-sectional images, one from each of the two |i-CT data sets, are shown in 
FIGs. 8(a) and 8(b). The two CT image sets were processed to generate BVF maps. Because of 
the small voxel size, the images have bimodal intensity distributions and could thus be 
segmented by selecting a threshold (1x77/) at the midpoint of the two modes representing the 
mean bone and the mean marrow intensities. Subsequently, the largest 26-connected 
component (Saha et a/., CVGIP: Image Understanding 63:418-429 (1996)) within the region 
thresholded for bone was computed. The intensities at the two modes were selected as the mean 
intensities for bone and marrow (denoted by and ixm respectively). At each point />, the BVF 
value was computed through the following equation: 

- 25 - P-2944 

520631 I 



Attorney Docket 22253-73413 



0, 



if/(/>) < 



if Mm ^f(P)<MB^ 



(Equation 14) 



Mb -Mm 

1, 



otherwise, 



where / denotes the image intensity function. 3D renditions of the B VF images generated from 
the two 3D |i-CT data sets of trabecular bone are presented in FIGs. 8(c) and 8(d), The surface- 
rendered 3D images were obtained using 3DVIEWNDC supported shell rendering (Udupa et aL, 
1993). 

[0104] For this and the following experiments, the bone sample displayed in FIG. 8(c) as 
sample A, and the cross section in FIG. 8(a) as slice A, are referred to herein as "sample A;" and 
those displayed in FIGs. 8(b) and 8(d), as sample B and slice B, respectively are referred to 
herein as "sample B." 

[0105] For BVF Object Generation Experiment 1 , the BVF images of the two bone 
specimens were resampled at voxel sizes corresponding to integral multiples of the parent 
resolution, i.e., 44 jam, 66 jam, ... 176 |im isotropic voxel size. This strategy was chosen for 
generating trabecular bone images at lower resolutions, as opposed to re-sampling the raw )x-CT 
data sets followed by computing the BVF images as the latter strategy would cause a significant 
loss in structural detail in the segmented images using the simple threshold-based method as 
used here. Segmentation may be improved by using more sophisticated methods Udupa et al, 
1996; Saha et a/., 2000); however, segmentation is not a focus of the present invention. 
However, for Experiments 2 and 3 of this Example, the respective operations were performed 
on the raw //-CT images first, and then the BVF computation method was applied. 

Experiment 1 - Performance of the Invention in Various Resolution Regimes. 
[0106] This first experiment was designed to assess the performance of the FDT method in 
various resolution regimes ranging from ex vivo |x-CT to in vivo |x-MRI data. Resolution was 
varied isotropically from 22 |Lim to 176 jun, in steps of 22 |imi, yielding eight different resolution 
values. FIGs. 9(a) and 9(b) display slices A and B from the BVF images of samples A and B at 
22 ^m isotropic voxel size. Matching slices from the BVF images at 88 jmi and I16\i m 
isotropic voxel size are displayed in FIGs. 9(c)-9(f). 
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[0107] Thickness values computed for the two different |a-CT image sets at different 
resolutions were plotted in FIG. 10 before and after the resolution dependent corrections 
(Equations 10-12). In spite of the wide range of resolutions examined, the corrected thickness 
values were shown to be stable as shown, with a maximum variation of 7.6% for sample A, and 
6.7% for sample B. 

[0108] A significant error was noted, however, in the computed thickness without correction 
(up to 45.3% for sample A and 49.8% for sample B). The error is negative, i.e., without the 
correction the algorithm underestimates true thickness values as predicted. However, there is a 
small positive error in the corrected thickness values, which increases approximately linearly 
with voxel size. This overestimation may result fi-om the artifactual filling of small lacunae (for 
example, see the apparent cavity shown by an arrow in FIGs. 9(a) and 9(c), which is filled in for 
FIG. 9(e). Likewise, merging of closely spaced trabeculae (for example, the region highlighted 
by an arrow in FIGs. 9(b), 9(d) and 9(f)) increased in apparent thickness at such locations. 
Equation 1 1 assumes structural anisotropy in the transverse plane. However, in reality there 
may be small deviations from transverse isotropy, meaning that the correction factor of 
Equation 1 1 is only approximately valid. 

[0109] It should be noted from the data plotted in FIG. 10 that the present method is 
remarkably reliable at voxel sizes (e.g., 176 jum), exceeding trabecular thickness (approximately 
124 jum and 110 //m). FIG. 10 also shows the computed thickness values using the HR method 
(Hildebrand et aL, 1997). Although at 22 //m, the estimated thickness values derived from the 
two methods are close, the disparity increases at lower resolution. It is noted that the HR 
method completely fails at voxel sizes greater than mean trabecular bone thickness. Moreover, 
the smallest thickness value at any location in any image computed by HR method is the voxel 
size itself Such behavior appears to be true for any threshold-based method. 
[0110] Further, even at 22 jum resolution, the change in computed thickness was 24% and 
32% for samples A and B, respectively, as the BVF threshold was varied from 25% to 75%. 

Experiment 2 - Rotation Dependence 
[0111] This experiment demonstrates the performance of the present method under object 
rotation. Thickness values were computed for the trabecular bone networks after the images 
were rotated in-plane, anticlockwise in steps of 5*^ over a range of 0® to 360*^ around the z-axis 

(axis perpendicular to image plane). Ge {15,10, ... 355} degrees were applied to each raw //- 
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CT image, and the rotated images were digitized in the original cubic grid using hnear 
interpolation. Images of a slice from sample B for three different rotation angles are illustrated 
in FIG. 11, and the results of this experiment are plotted in FIG. 12. 

[0112] BVF images were computed using Equation 14. Notably, the variations in measured 
thickness do not exceed 3.3% and 4.5%, for samples A and B, respectively. Both rotation- 
thickness plots show a periodicity at integral multiples of 45®. This observation is consistent 
with the behavior of digital cubic grids in which effects of rotation are largest around 45® and its 
integral multiples. 

Experiment 3 - Noise Dependence. 
[0113] This experiment was designed to examine the sensitivity of the present method to 
noise. FIG. 13 displays slice A for three different SNR values, while plots of apparent thickness 
versus SNR are illustrated in FIG. 14. Zero mean, multiplicative random noise corresponding to 
SNR e {1, 2, ... 20} was added to the //-CT images of samples A and B. The SNR value was 
calculated as the ratio between the mean and standard deviation of intensities over a sufficiently 
large marrow region. As in the previous Example, BVF images were computed using Equation 
14. 

[0114] Notably, at very low SNR values (<5), computed thickness values were substantially 
underestimated (almost 50%). However, above SNR « 5 the values stabilized, asymptotically 
approaching their noise-free values. For SNR > 5, variations in estimated thickness value for 
the two samples were less than 3.5% and 2.7%, respectively. The likely cause for the method to 
break down at very low SNR is the creation of tiny cavities and lacunae in the support of 
trabeculae caused by noise, thus increasing apparent bone skeletal area, resulting in an 
artifactual reduction in computed thickness. 

Experiment 4 - Reproducibility. 
[0115] This experiment was designed to examine to what extent spatially close slice 
locations in the same specimen return similar thickness values. The hypothesis to test was that 
two similar locations within a volume of trabecular bone volume would yield similar thickness 
distributions (since, the trabecular architecture varies continuously along the slice direction). In 
contrast to the preceding examples, images were used only from individual slices. Specifically, 
images from two slice locations separated by 330 pim (15 slices) were selected from each of the 
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BVF images of samples A and B. The rationale underlying this experiment is the notion that the 
structxiral characteristics foimd within the bone of a particular donor had to be more similar than 
those in bone from different donors. 

[0116] As shown in FIG. 15, two pairs of FDT-based trabecular bone thickness distributions 
are displayed without resolution dependent corrections, each for the pair of slices from one of 
the two samples. 

[0117] Notably, from these plots that thickness distributions for the slices from the same 
sample are quite similar while differing substantially for the slices originating from different 
samples. 

Example 4 - Reproducibility of Method In Vivo 

[0118] A first experiment was aimed at evaluating the reproducibility of the proposed 
trabecular bone thickness computation method in human subjects using in vivo //-MRI data. 
Therefore, high-resolution images of the distal radius in four human volunteers were examined. 
Each individual had been previously scanned three times with repositioning and reset apparatus- 
setup as part of a project designed to evaluate the serial reproducibility of trabecular bone 
structural parameters. Image data sets consisted of 28 shoes of 137x137 fxrvL' pixel size having 
350 fixa thickness, obtained with the 3D FLASE partial flip-angle spin-echo sequence (Ma et 
a/., Magnet. Reson. Med 35:903-910 (1996); Song et al. Magnet. Reson. Med. 41 :947-953 
(1999)). For each image, an actual volume of interest (VOI) for which thickness was computed 
was generated by manually tracing boimdaries on each slice using an IDL-based graphical 
interface running on a Pentium III PC under MS Window OS. The data encompassed by the 
VOI were preprocessed by deshading and noise reduction using a histogram deconvolution 
method (Hwang et al., 1999) to produce a set of BVF images. SNR values ranged between 12 
to 14. 

[0119] An unprocessed //-MR image is shown in FIG. 16(a), and a zoomed region of 
interest fi-om the same slice is displayed in FIG. 16(b). The boundary of the VOI, used for 
thickness computation, for the particular slice is also shown FIG. 16(a). A matching slice image 
fi-om the BVF data is shown in FIG. 9(c), along with the computed FDT map (FIG. 16(d)). 
[0120] The results of these experiments are important because they determine the power of 
the method in longitudinal studies, designed, for example, to evaluate the efficiency of 
treatment. The reproducibility data yielded trabecular bone thickness means ranging from 93 to 
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124 fim. This analysis resulted an average coefficient of variation of 3.0% and an intra-class 
correlation coefficient of 0.96. Such precision allows monitoring of progression or regression of 
disease in serial drug intervention studies as discussed below. 

[0121] A second experiment was directed toward to evaluating the performance of a method 
of the present invention in a longitudinal study designed to evaluate the effect of supra- 
physiological levels of corticosteroids on trabecular bone. Detecting the effect of steroid 
treatment was selected because it is a well known cause of thinning of trabecular bone. 
[0122] For this purpose, data fi-om a study conducted on rabbits which had received 
dexamethasone (a synthetic analog of Cortisol) as a means to induce bone atrophy (Takahashi et 
al, Proc. Natl Acad, Scl USA 19:19 (2002)), were re-examined. In brief, the protocol involved 
eleven New Zealand white rabbits that either received dexamethasone (0.4 mg/kg/day, N = 6), 
or were sham operated (N = 5). Images were obtained at three time points, first at baseline and 
subsequently at four (4) and eight (8) weeks after the beginning of the corticosteroid-treatment. 
MR images of the distal femur epiphysis were acquired with the FLASE pulse sequence used in 
the previously described human subject study (Example 7), affording 28 contiguous slices of 
97x97x300 |im^ voxel size. The visually matched images of |i-MRI data sets (SNR values 
between 9 and 1 1), as illustrated at FIG. 17, were subsequently processed to yield BVF maps 
using a method analogous to the one used for human radius images. 
[0123] The hypothesis tested was that the corticosteroid-treated rabbit group would lose 
bone due to thinning of trabeculae relative to the sham-operated (control) animals. The results 
of FDT-derived thickness analysis of image data collected previously (Takahashi et al.y 2002) 
are plotted in FIG. 18. The plots show the mean average thickness values (after resolution- 
dependent correction) for each group at each of the three time points - baseline, 4-weeks, and 8- 
weeks fi-om onset of treatment. The data confirm the short-term effect of corticosteroid 
treatment on trabecular bone, indicating relative decreases in trabecular bone thickness of 2.8% 
at 4-weeks and 8.4% at 8-weeks. By contrast, however, trabecular thickness increased in the 
sham-operated group by 0.55% at 4-weeks and 5.3% at 8-weeks. The differences between 
treated and controlled groups were statistically significant (p 0.05 at 4 weeks and p = 0.02 at 8 
weeks). 

[0124] Although bone loss in humans may be slower than in animal models, trabecular 
thinning subsequent to glucocorticoid treatment can be quite fast. As the data above show, bone 
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loss and the concomitant reduction in trabecular thickness can be extraordinarily rapid. Studies 
on the basis of histomorphometry on iliac bone biopsies in patients on short-term glucocorticoid 
treatment exhibited a reduction in trabecular thickness of about 20% relative to age-matched 
women with postmenopausal osteoporosis (Carbonare et al, J. Bone Miner. Res. 16:97-103 
(2001)). Similarly, short-term administration of parathyroid hormone was foxmd to substantially 
increase trabecular thickness (Bradbeer et al, Clin, Endocrinol (Oxf) 37:282-289 (1992)). 
Thus, the effects detected in the present studies are of the same order as those presented above 
in a rabbit model of steroid-induced osteopenia, supporting the conclusion that the present 
methods provide an improved capabiHty to follow patients during treatment. 
[0125] These observations were confirmed by the mean thickness values. Thus, as shown 
by the disclosed experiments, the thickness value computed using the methods of the present 
invention is an intrinsic property of trabecular bone at a particular skeletal location in a subject, 
although these particular values may be differ for corresponding locations in different subjects. 
[0126] The disclosures of each patent, patent application and publication cited or described 
in this document are hereby incorporated herein by reference, in their entirety. 
[0127] While the foregoing specification has been described with regard to certain preferred 
embodiments, and many details have been set forth for the purpose of illustration, it will be 
apparent to those skilled in the art without departing firom the spirit and scope of the invention, 
that the invention may be subject to various modifications and additional embodiments, and that 
certain of the details described herein can be varied considerably without departing fi*om the 
basic principles of the invention. Such modifications and additional embodiments are also 
intended to fall within the scope of the appended claims. 
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